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We study a random heteropolymer model with Langevin dynamics, in the supersymmetric formu- 
lation. Employing a procedure similar to one that has been used in static calculations, we construct 
an ensemble in which the affinity of the system for a native state is controlled by a "selection tem- 
perature" To. In the limit of high To, the model reduces to a random heteropolymer, while for 
To ^ the system is forced into the native state. Within the Gaussian variational approach that 
we employed previously for the random heteropolymer, we explore the phases of the system for high 
and low To. For high To, the system exhibits a (dynamical) spin glass phase, like that found for the 
random heteropolymer, below a temperature Tg. For low To, we find an ordered phase, characterized 
by a nonzero overlap with the native state, below a temperature T„ oc 1/To > Tg. However, the 
random-globule phase remains locally stable below T„, down to the dynamical glass transition at Tg. 
Thus, in this model, folding is rapid for temperatures between Tg and T„, but below Tg the system 
can get trapped in conformations uncorrelated with the native state. At a lower temperature, the 
ordered phase can also undergo a dynamical glass transition, splitting into substates separated by 
large barriers. 
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I. INTRODUCTION 



The protein folding process is relevant for all aspects of 
life: once read off from the RNA chain, proteins perform 
a variety of functions, from mechanical work to attacking 
viruses. [1] The key factor which determines the function 
of a protein molecule is its 3D structure, which, in turn, 
is determined by the sequence of amino acids forming the 
protein chain. [2-5] Furthermore, a protein that has been 
denatured (by stretching it for example) finds its native 
state relatively quickly. Protein folding has attracted an 
enormous amount of scientific attention, but still there is 
no generic understanding of this process. Nevertheless, 
one thing is clear: a proteins generally has a potential 
energy surface which results in a stable free energy min- 
imum, corresponding to the native state [3]. 

Random heteropolymer models (RHP) have been used 
extensively as candidate systems which might help us 
to understand the generic features of the potential en- 
ergy surfaces of proteins and their connection with ther- 
modynamic [6-13] and dynamical [14-20] properties. 
The RHP model is characterized by quenched random 
monomer- monomer interactions, meant to mimic the va- 
riety of interactions between amino-acids in random se- 
quences. It turns out that the potential energy surface 
of the RHP is quite similar to that of a particular class 
of spin glasses [21]: Its complex form, with exponentially 
large numbers of local minima and saddle points, con- 



strains the motion of the system drastically, and it can- 
not explore its full configuration space and reach Gibbs 
equilibrium. In a previous paper ( [22], henceforth re- 
ferred to as paper I), we demonstrated, in mean field the- 
ory, the existence of a sharp transition to a "dynamical 
glassy state" in which the equilibration time diverges and 
the dynamics exhibit aging. (The potential importance 
of spin glass physics to proteins was first discussed in 
Ref. [23]). Obviously, the random heteropolymer model 
does not describe a protein with a native state, but it 
alerts us to the need to examine possible glassiness in 
models for protein dynamics. 

Why are real proteins not glassy? Evidently, nature 
has tuned amino acid sequences to avoid glassy behav- 
ior. To understand how such tuning might be done, it 
is worthwhile to study models which contain competi- 
tion between glassiness and a tendency to form a native 
state, by choosing interactions which are not completely 
random. Several studies along the lines of this suggestion 
have been made in statics (using the replica treatment, 
see, e.g., Ref. [23]). The tendency toward a particular 
state can be built in by choosing sequences from a dis- 
tribution correlated with the native sequence [2,24-26]. 
A dynamical treatment of similar models is highly desir- 
able, not only to help gain insight into results obtained 
in replica approaches, but also because knowledge of the 
correct thermodynamics alone may not be sufficient: it is 
known that in related (mean field models) static and dy- 
namic phase diagrams can be different. Thus (at least on 
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sufRciently short time scales) only a dynamical approach 
can describe the measurable properties of the system. In 
this paper we undertake such a study. 

Wc extend the RHP model studied in [6,7] to include 
the existence of a native state: the original random 
monomer-monomer interactions are biased so as to favor 
the native state conformation. The problem is formu- 
lated as a Langevin model. To the best of our knowledge, 
there is so far neither a static nor a dynamic treatment 
available for a model of this sort: Static studies have been 
based on random monomer sequences, i.e., using only N 
random parameters, see Rcfs. [2,24-26], rather than the 
N{N - l)/2 in the RHP model. 

Admittedly, the model does not describe a realistic 
protein (e.g., it does not give rise to secondary structure 
such as a-heliccs or /3-sheets). However, it docs con- 
tain important generic features: the polymeric structure 
and the mixture of attractive and repulsive interactions. 
Together, these fcatmes lead to frustration in the struc- 
tural dynamics. In our view, ours is the simplest such 
model that includes competition between glassy and na- 
tive states. As we will see, it teaches us that one can not 
get rid of glassiness so easily. 

As in paper I, we simplify the model further by omit- 
ting three-body interactions in the polymer. (A re- 
view describing how to include three-body terms is given 
in [9].) The price we have to pay for this simplification 
is that wc have to introduce a somewhat arbitrary con- 
fining potential, which we take to have a quadratic form. 
We adjust its strength so that the radius of gyration Rg 
of a polymer of size N scales like N~^/'^, where d is the 
dimensionality of the system. In this way we attempt to 
describe a globular state. Of course, we can not describe 
the ^-point transition in such a model, but here we are 
only interested in transitions between different globular 
states. 

Our formal starting point is the Martin-Siggia-Rosc 
generating functional for the Langevin dynamics of the 
model [27-30], written, for convenience and compactness, 
in its supersymmetric form [31]. To derive equations 
of motion for correlation and response functions wc use 
a variational ansatz with a quadratic action. This ap- 
proach has been used to study the problem of a manifold 
in a random potential, in both statics [32,33] and dynam- 
ics [34,35]. 

In paper I we showed that the RHP model exhibited 
broken ergodicity (formally, a spontaneous supersym- 
metry breaking) in a low-temperature dynamical glassy 
phase. In the present study, with interactions biased 
in favor of a native state to a controlled degree, wc 
find, in addition, a well-folded phase, if the bias is 
strong enough. It can coexist with either the disordered 
(random-globule) state or the frozen-globule glass phase, 
depending on the temperature. Furthermore, we find 
that at low temperature the native phase can itself un- 



dergo a dynamical freezing into a different glassy phase. 
In this phase the conformation of the protein is always 
highly correlated with the native state, but cooperative 
kinetic constraints still lead to a divergent equilibration 
time, as for the frozen-globule state. 



II. THE MODEL 

The model is defined as follows. The Langevin dynam- 
ics is assumed to be governed by a Hamiltonian H[x], 

dx{s,t)/dt = -5H[x]/6x{s,t) + ■n{s,t). (1) 

Here x(s,t) is the position of monomer s at time t. The 
monomers are numbered continuously from s = to 
8 = N. r]{s, t) is Gaussian noise 

{r,{s, t)r,{s', t'))T = 2TS{s - s')5it - t'), (2) 

resulting from coupling to a heat bath at temperature T. 

The Hamiltonian H[x] contains a deterministic part 
Hf)[x,fj] and a random part H[x,{B}]. Ho[x,ij] is de- 
fined as 

Ho[x,tj] = - ds{[dx{s,t)/ds]''+tix{s,tf}. (3) 
^ Jo 

It describes the elastic properties of the chain and a con- 
finement potential which fixes the density of the pro- 
tein. The radius of gyration Rg ~ so, in order 
that the protein is compact, i.e., Rg ^ N'^^'^, we require 
jjL ~ A/'"*/'^. Thus, since we are interested in very long 
proteins (to obtain the thermodynamic limit) we need to 
solve the model for /i close to zero. 

The random part H[x, {B}] describes the quenched 
random interactions between monomers, 

1 

H[x,{B}] = - dsds'B,,,V{x{s,t)-x{s',t)). (4) 
^ Jo 

We take Bgs' Gaussian, with variance B^. The 
quenched average over Bss' is performed as ((.))s = 
/ns>s' dBss'{-)P{{B}). V{Ax) is a short-range poten- 
tial, and, for simplicity, we take it to have a Gaussian 
form, as in Ref. [16], 

V{Ax)=(-^] e-(^^)'/2._ 

\ 2770" / 

d is the dimensionality of the system and y/a the range 
of the potential. Large (small) a corresponds to a 
long (short) range potential. In particular, for cr ^ 0, 
V{Ax) S{Ax), and we recover the potential used 
in [6,7,18]. Here and in the following Aa; refers to a 
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monomer-monomer distance: Aa; = x{s,t) — x{s',t) for 
a pair of monomers s, s' . 

We use reasoning similar to that employed in statics 
to define P{{B}) (see Refs. [2,24-26]), adapting it to the 
random-bond model: 



P{{B}) a e"^^'^''''^^"^ J'isds'Bl,/2B^ 



(6) 



To is called the selection temperature, and Xq{s) is some 
arbitrary native state conformation. Thus the symmetric 
bond distribution of the RHP model is distorted so as to 
give bigger weight to B^s' 's which are attractive between 
monomers which lie close to each other in the configura- 
tion Xo{s). Explicitly, the properly normalized P{{B}) 
is given by 

P{{B}) = (27rS2)"^(^~^)/'*e"'^o"^^^'* / dsds'v(xo(s)~xo{s')f 

^g-/3o/2 J dsds'B,,,V(xo(s)-xo{s'))-l/2 J dsds' B^^, /2B^ 

from which we see that the distribution of Bss' is peaked 
around S™?^ = - (3oB^V {xq{s) - xo{s')). Thus, if 
monomers s and s' are close in the native state {V{xo{s) — 
Xo(s')) 7^ 0)> their coupling constant Bg^' is pulled down, 
as in a Go model [36,37]. For Tq — > oo we recover the 
RHP model. For Tq 0, P{{B}) picks a specific set 
of Bss'- For this set, by construction, xo{s) is the deep- 
est minimum of H[x, {B}] given in Eq. (4). This is the 
mechanism that embeds the native state xo{s). 

This mechanism is somewhat arbitrary. However, 
the fact that the strength of embedding of the native 
state is controlled by the single parameter To facilitates 
the study of transitions between random and native-like 
states (and, as we will show, of possible coexistence of 
such phases). 

So far, the configuration a;o(s) is arbitrary. Thus xo(s) 
has to be considered a quenched random function, to be 
averaged over just like Bss' in order to obtain generic 
results. We will carry this average out later. 

All our results are obtained in the thermodynamic 
limit, where the length N of the heteropolymer chain 
goes to infinity. Also, for simplicity, we join the polymer 
ends to form a ring. This neglect of end effects is valid 
for a long chain. 



III. MAPPING TO THE FIELD THEORY 



To solve the model we map the Langevin dynamics 
onto a supersymmetric (SUSY) field theory. Using the 
standard Martin-Siggia-Rose formalism [27-30] and su- 
persymmetric (SUSY) notation [19,20,31,38], the dynam- 
ical average of any observable, for fixed {B}, can be cal- 
culated as (see, e.g.. Paper I for details). 



{0[^])t = J £'$0[$]e-^[*l, 
5[$] = 5i[$]+5[$,aro,{B}], 



where 



(8) 
(9) 

(10) 



S'i[$] = l/2 j dsdlds'd2^s,l)Kt'2'Hs'2), 

S[^,xo,{B}] = 1/2 j dldsds' x 

xBs,s'Vms,l)-^s',l)), (11) 



and 



i^ri = Si2Sss'Kt , Kt=T[ii- {d/dsf] - 

^2 d 



Df)=2T ^ 



26*1- 



ddidOi dOidti dti 

The $(s, 1) denotes a superfield 

^{s,l)=x{s,h) + e^ri{s,ti) + 

+f]{s,ti)ei+0i0ix{s,h) 



(12) 
(13) 



(14) 



containing the physical coordinate x{s,t), the MSR aux- 
iliary field x(s,t), ghost fields rj^s^t) and ??(s,t) that en- 
force the normalization of the distribution, and Grass- 
mann variables 6 and 9. We use the notation 1 = 
(6'i,^i,ii), likewise J dl = J d9id9idti. 

Of course, the solution can be obtained without the 
aid of the supersymmetric formalism, but we find it con- 
veniently compact. 

As noticed by De Dominicis [28] the expression in 
Eq.(8) is already normahzed, so the average over the 
quenched disorder Bs^s' can be done directly on (8): 

{{0[^])t)b = J T)$C>[$]e-(^iI*l+^^[*'^°l\ (15) 
where exp(-S'2[$, xo]) = {exp{-S[^,xo,{B}]))B, and 
S2[^, xo] = J dsds' J dlVms, 1) - $(s', 1)) 



xV{xo{s) - xo{s')). 



(16) 



Thus, the native state xois) enters the action in the 
second term of Eq. (16). Note that there is no term 
(3qV{xo{s) — xo{s'j)'^, since it gets cancelled by a similar 
normalization factor for P({B}) in Eq. (7). It is useful 
to rewrite Eq. (16) as 

52 = -^ J d'^xd'^ydl d2A'^^\x,y)A[^^{x,y) 



j d^x d^y dl a(^) {x, y)A['^ {x, y) (17) 
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with the notation A[{\x,y) = Jd.sf{^{s,l) - 

x)fms, 2)-y), A[i^ {x, y) = J dsfi^s, 1) - x)f{xo{s) - 
y); f & In the long-chain Umit, as discussed 

in Paper I (and references therein), one obtains a self- 
consistent field theoretic formulation, with 5*2 simplified 
to, 

S2[<^,xo] = ^J <fx<fydld2 [{A^^^\x,y)){A^^{x,y))- 
-A'^^^\x,y){A^^{x,y)) - {A^^^\x,y))A^^{x,y) 



-A'^,\x,y){AY^{x,y)) - {A\l> {x,y))A\y {x,y)\ . (18) 

All averages of the type {A'-^'^^) have to be calculated 
self-consistently with S'[$] = + S'2[$]. (We have 

abbreviated the double average {{■)t)b simply by (.).) 
In the limit N ^ oo Eqs. (15) and (18) provide an exact 
description of the dynamics for an arbitrary native state 
xo{s). 



d'^xd'^ydl 



,(-5) 



{AZ\x,y)){A^^{x,y)) 



{V)i 



IV. AVERAGE OVER NATIVE STATE 
CONFORMATIONS 



It is impossible to solve the model for a general na- 
tive state configuration xo{s). We therefore consider a 
distribution of native states and perform the average 



< 0[^,xq\ > = j Dxo < 0[^,xo] > 



->S'o[a;o] 



(19) 



where 5*0 [a;o] weights each native state conformation in 
the ensemble as 



! J dsxo 



with 



So[xo] = 1/2 / dsxo{s)K^^ xois') 



(20) 



(21) 



The parameter fiQ fixes a size of the globule in this en- 
semble, 



2Vmo 



(22) 



Since the polymer ends are joined, there is translational 
invariance along the coordinate s and {xo{s)^) does not 
depend on s. Thus, with this procedure, the dynamical 
generating functional for the problem is calculated as 

^-F^yn = J £,a.Q£,$e-(So[xo] + Si[4.]+S2[*,Xo])_ (23) 



There is some formal similarity between the dynamical 
functional -Fdj/n a-nd the static replica partition function. 
The integration over Dxq enters in the same way as the 
extra replica in the static formalism. 



V. CORRELATION FUNCTIONS 

The SUSY correlation functions 

G'l^' = ($(s,l)$(s',2)) (24) 

GTo ^ ms,l)xo{s')) (25) 

^ (.To(.s).To(.s')> (26) 

contain all the information we are interested in. 

Gfl encodes 16 correlation functions, out of which only 
two, correlation and response function, are independent 

and nonzero: 



G 
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X [02R{s, ti;s', t2) - 9iR{s', t2; s, h)], (27) 



with 



C{s,t;s',t') = {x{s,t)x{s',t')), 



(28) 



R{s,t;s',t') ^ {x{s,t)x{s',t')) = ^1^. (29) 

The field h{s',t') entering the description of response 
function is an arbitrary external field that couples to 
x{s', t'). The fact that only two correlation functions sur- 
vive is related to Ward identities originating from SUSY 
invariance of the original action S. 

The supersymmetry of the theory is associated with 
equilibrium. One of the Ward identities resulting 
from SUSY is the fluctuation-dissipation theorem (FDT) 
which relates correlation and response functions. In the 
present case, the glassy state manifests itself as a sponta- 
neous breaking of supersymmetry, leading to a modified 
FDT, as in previous treatments of other models [31,34]. 

Gfo describes the overlap with the native state. Due 
to Ward identities, only a single correlation function sur- 
vives (see Appendix A for details): 

G\i = {x{s,t)xo{s')) = <P{s,tv,s'). (30) 

Similarly, the native state ensemble is described by 

G'ot; = {xois)xo{s'))=r{s;s'). (31) 

Gfl alone is suSicient to describe the RHP model. Here 
we need the two extra functions Gf g and Ggo • 

Also, in what follows, we exploit the translational 
invariance along the s coordinate and define Fourier 
transforms of all correlation functions: X{s,s') = 
J ^^iKs-s')x^ where X = C,R, </., T. 
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VI. EQUATIONS OF MOTION 

To solve the model we proceed by making a Gaussian 
variational ansatz (GVA), assuming that the fields $ are 
described by the approximate action 

Svar = \ j dldsd2d8'^{s,l){G-^)li^{s',2) + 

+ j dldsds'^s,l){G-^){ixo{s') + 

+i j dsds'xo{s){G-^)lixo{s'). (32) 

Technically, this implies the following approximation for 



where 



The stationarity condition 



SG 



(33) 

(34) 
(35) 

(36) 
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translates into the equation of motion for Green's func- 
tion Gfl (see Eqs. 43-46). We have derived iden- 
tical equations of motion by using the approach of 
Ref. [35], where standard field theoretic identities (e.g. 

= ) are used. It can be shown that for 
quadratic S^ar the two procedures give the same result. 
We omit this analysis here to save space. 

In a corresponding equilibrium problem, the stationar- 
ity condition is also an extremum condition and provides 
a bound on the free energy. Here, since F^yn contains in- 
tegrations over complex fields and Grassmann variables, 
the GVA docs not give a bound on Fdyn- Nevertheless, it 
is the first step in a systematic approximation procedure, 
as outlined in Appendix B. 

The GVA has been applied to the problem of a mani- 
fold in a random potential, in both statics [32,33] and dy- 
namics [34,35]. The method is exact when the dimension- 
ality of the manifold is infinite but is only approximate 
for finite dimensionality. Nevertheless, even for rather 
low dimensionality it has been shown to be a very good 
approximation in the random-manifold problem, where 
it has been checked numerically [35]. We have shown in 
Paper I that the present model is closely related to the 
random-manifold problem. Thus, we hope that the GVA 
will also be reasonable here, although we have not strictly 
checked its validity. 



Using (33), (32) and (9) gives the following expression 
for Fdyn- 



- I dsds'K^^ GgS + 



^ f dsds' dld2Kl( G'^^ - ^TrlnG 
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dVydld2{A^^^\x, y)){A^^{x, y)) 

j d''xd''ydld2{AZ\x,y)){A'^^{x,y)), (37) 



where all averages are to be calculated using S^ar (see 
Eq. 32). Performing averages, the fourth and fifth term 
on the right hand side of (37) become 



F, 



(4) 



dld2dsds'V 



F^ll = -^ I dld2dsds'V 



{BI2 + B{2)/2_ 
\bI^ + b(^)I2 



(38) 
(39) 



where 



B{2 - 1) - 2)]') = Gfi + - 2Gr2, (40) 

Bi, = ms, 1) - xo{s)f) = Gfi + Gil - 2Gfo, (41) 



and 



V{z) = -^{z + a)-'''^ , B^ = ^-{4n)-^/^. (42) 



Performing the variational ansatz (i.e., evaluating 
Eq. 36) results in the following equations of motion: 



Till + e) - D\ 



(2) 



^12 = ^12 + 2 



j d3V'{B,s) 



x(G^2 - GI2) + 2poV'iB,o){G'o2 - 



(43) 



Till + e) - Dp^] G^o = 2 y" d2V'{Bi2) X 
(G^o - G\^) + 2/3oV'(Bio)(Ggo - G^o), (44) 



(mo + fc2)Ggi = 2/3o j d2V'(B2o)(G^i - Ggi), (45) 
(mo + A;2)GSo = l + 2Poj dlV\B,o){Glo - G^o), (46) 



and after disentangling the SUSY notation one gets (see 
Paper I for related details) 



5 



[T{ti + k'') + d/dt]Ck{t,t') = 



and, in the aging regime, 



2TRk{t', t) + 2 I dt"V [B{t, t")] Rk{t', t") 
Jo 

+4 / dfV'[Bit,r)]rit,t")[Ck{t,t')-Ck{t",t' 
Jo 



-20oV'[AmCkit,t')-Mt% 



(47) 



[T{n + + d/dt]Rk{t, t') = 5{t - t') + 

+4 /* dt"V" [Bit, t")] r{t, t") [Rkit, t') - Rkit",t')] 
Jo 

-2poV'[A{t)]Rk{t,t'), (48) 

[T{f, + k'') + d/dt]Mt) = 

4 /* dt"V" [B{t, t")] r{t, t") [<t>k{t) - (t>k{t")] 
Jo 

+2poV[A{t)]{Tk-<i>k{t)), (49) 



(/xo + fc^)0fe(i) =2/3o / dt"V[A{t")]Ru{t,t"), (50) 
Jo 



(mo + k')Tk = 1, 



(51) 



where J5(i, t') and A{t) arc defined as B(t, t') = {{x{s, t) — 
x{s, t') f) = C{s, t; s, t) + C{s, t'; s, t') - 2C(s, t; s, f) and 
A{t) = {{x{s,t) - xo(s))2) = C{s,t;s,t) - 2(f>{s,t;s) + 
r(,s, s). Note that due to translational invariancc with re- 
spect to s both B{t,t') and A{t) are s-independent. The 
equations of motion for Ck{t,t') and Rk{t,t') are almost 
identical to those for the pure RHP model. Coupling to 
the native state enters through the terms proportional to 
I3q. Again, for large selection temperature /3o — > and 
one recovers the RHP model. 



lim Ck{t, Xt) = qkCkiX), 

t^OO 

lim Rk{t,Xt) = 7i?fe(A). 



(54) 
(55) 



The validity of these assumptions could be checked nu- 
merically. Since this has been done for equations of simi- 
lar type elsewhere [35], we omit it in the present analysis. 

Second, it is well known that asymptotic solutions of 
such equations can be characterized by a few order pa- 
rameters [34,35,39^1]. They are defined as 



qk = lim Ck{t,t), 

t— ^oo 

qk = lim Cfc(T), 



go,fc = lim^QkCkiX), 

t— ^oo 



(56) 
(57) 
(58) 
(59) 



The following fe-integrated quantities will also be useful: 
dk 



(60) 



/dk 
—qk= lim lim {x{s,t)x{s,t + T)}, (61) 
ZTT T—*oo t— ^oo 

/dk 
— (7o.fc = lim lim {x{s,t)x{s,Xt)), (62) 

/dk 
—ipk= lim {x{s,t)xQ{s)). (63) 
iTT t— >(X) 



q measures the size of the globule, q measures the per- 
sistent correlation in the TTI regime, go the asymptotic 
correlation in the aging regime, and the overlap with 
native state. Also, it is useful to define 



VII. EXTRACTING ORDER PARAMETERS 



The equations of motion are coupled integro- 
differential equations with initial conditions given by 
Cfe(0, 0), 0(0), and (we use Ito's convention) R{t + e, t) 
1 as e ^ 0. To solve the equations analytically we have 
to consider several assumptions (which can be checked 
by numerical solution). 

First, we make the (rather strong) standard assump- 
tions from aging theory for spin glasses about the asymp- 
totic behavior of the solutions: In the regime of time 
translational invariance (TTI), 



lim Ck{t + T,t) = Ck{T), 
lim Rk{t + T,t) = Rk{T), 

t^OO 



(52) 
(53) 



b = 2{q-q), ho = 2{q-qo), 



(64) 



a = lim {[x{s, t) - Xo{s)f) =q-2^+ ——. (65) 



Third, we assume that the generalized fluctuation dis- 
sipation theorem is valid in the form 



^ a; dCk{X) 



(66) 



X could in principle depend on k and Ck- However, re- 
lated models have been studied in detail and they exhibit 
one step replica symmetry breaking with a fc-indcpcndent 
X. This one step replica symmetry breaking ansatz in our 
dynamical study translates exactly to Eq. (66). 
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VIII. RELATING ORDER PARAMETERS 



For t = t' and t — > oo Eq. (47) gives 
T(m + = T + |V'(6)(1 - x){qk - qk) 

+ ^V'{bo)x{qk - qo,k) - 2/3oV'(a)(9fc - (pk). (67) 



With t = t' + T and t' ^ oo and then r ^ oo Eq. (47) 
becomes 

T{i^ + k')qk = ^{V'{b) - xV'{bo)){qk - qk) 

+ ^V'{bo)x{qk - qo,k) - 2poV'{a){qk - <l>k)- (68) 

Eq. (47) in the aging regime t' = Xt, first for f — > oo and 
then A ^ 0, gives 

Till + k^)qo,k = |;V'(6o)(l - x){qk - qk) 

+ ^Vibo)x{qk - qo,k) - 2l3oV{a){qo,k - h)- (69) 



bo 



1 1 



x-1 1 



X X + 



. bo 
9 = TT 



V'ibo 



1 



1- ^ 



2 4T2/i3/2 4^\^i_i^^ 
x(2+W^)(l 



bo ^ V'jbo) ^ 1 



2 4T2/i3/2 4Vm^ (l + 



/MO ^^^2 



,/f /Mo 



(75) 



(76) 



(77) 



(78) 



Eqs. (49) and (50) result in two equations for <^fe, 

Tin + k^)<fik = ^V'{a)iTk - <Pk) (70) 
-to 
2 

(/zo + k'^)(Pk = 7^^'{a)3:{qk - qo,k) 
-t Jo 

+ 7^V'{a){qk - qk) (71) 

They are equivalent; one can chose to solve for the order 
parameters working with either (70) or (71). This seems 
a rather remarkable coincidence. We believe that it orig- 
inates from the SUSY invariance of the original action 
S. For example, a similar comment holds for equations 
(47) and (48) ; they arc equivalent in the TTI regime and 
one can derive one from the other. The 'conspiracy' of 
(49) and (50) not contradicting each other is very likely 
a similar phenomenon. Eq.(48) for A = 1 reduces to 

4V"(h) 

Rkim + k' + ^)^ -{qu - gfe)^-iif(l), (72) 



rp2 



where 



r(A) ^ / ^^Rk{X) 



and S is defined by 



^ = x^ (V'(6) - V'(6o)) . 



(73) 



(74) 



M = M+^V'(a), 



(79) 



and the combination of Eq. (75) and (72) gives 

Q = f{l)[T^ + b^V"{b)]. (80) 

Furthermore, the overlap <p with the native state is given 
by 



^ = 



1 1- ^ 

2jjio 1 4- 

V A 



(81) 



Solving these equations for the order parameters gives 



All overlap order parameters are positive. However, this 
result is not obvious and has to be obtained after some 
algebra. 

These equations have two kinds of solutions. In one 
kind, b = bo, so there is no glassiness (aging). For this 
kind of solution, the parameter x is irrelevant. We call 
such solutions "ergodic". (While it will turn out that 
some of them are not truly ergodic, in the sense of de- 
scribing states where the entire configuration space is vis- 
ited with Boltzmann probabilities, they violate ergodicity 
in a rather trivial way, like a ferromagnet below the Curie 
temperature. We could call them "non-glassy", but we 
prefer not to use a negative term.) 

For an ergodic solution, with b = bo, E = 0. Further- 
more, f(A) = 0, so Eqn. (80) is trivially satisfied. One 
then has to solve the four equations (75) and (77-79) for 
b, q, a and /i. 
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The stability of such a phase against glassiness can be 
determined using the analysis we presented in Paper I 
(see Fig. 1). There, we studied a model with no native- 
state bias in its interactions (Tq = oo) for finite ji. The 
boundary of the glassy state as a function of ^, has a form 
qualitatively like that in the p-spin glass as a function of 
field [42,43]. In the present model, the presence of the 
native state enters the calculation solely through the re- 
placement of /i by fi. Therefore, if a particular T and /}, 
fall in the glassy regime (the region below the full and 
dashed lines) in Fig. 1, the ergodic ansatz has to be given 
up. 

The instability can occur in two ways, according to 
whether /i is bigger or smaller than the critical value jlc- 
Above /},,;, the line separating glassy from ergodic regions 
is an Almeida-Thouless (AT) line; below it the stability 
condition 



T^ + b^V"{b) >0, 



(82) 



is violated. For jj, < jlc, there is no AT instability. The 
transition is like that for the completely random het- 
eropolymer. To find such a transition, we have to solve 
for a glassy phase, characterized in part by a value of 
the FDT- violation parameter x <1 and then find where 
in the parameter space x ^ 1. In the region where the 
a; < 1 solution exists, the associated ergodic phase is un- 
stable and is replaced by the glassy one. 

In a glassy phase, aging is present: f(l) ^ 0, so the 
quantity in brackets in Eq. (80) has to vanish, i.e., the 
AT condition has to be satisfied as an equality, rather 
than an inequality. This so-called marginal stability con- 
dition determines 6 as a function of temperature. In this 
case we have three more unknowns, S, bo and x, making 
a total of seven, and seven equations, (75-80), to solve 
for them. 

We look for ergodic solutions first in the next section, 
and we examine their stability. Then, in the following 
section, we study glassy solutions (within the 1-step ag- 
ing ansatz of section VII) and identify the regions in the 
parameter space where they hold. 



IX. ERGODIC PHASES 



For ergodic phases, Eqns. (75-79) reduce to 

b = bo = ^ (83) 

^ 2Vm ^4^o[l~fJ 



X 2 + 




1 ^'(7n^ 




j2 = ti+—V'{a) 



(85) 
(86) 



They can be solved numerically: given /i, ^o, T and Tq 
one can find /x, which in turn determines q, b = 6o (cquiv- 
alently q = go), and tp. However, it is possible to gain 
some analytic understanding in a few soluble limits. 

In this discussion we will concentrate on the limit of 

small /i. As we noted in paper I, if we want to confine 
N monomers within a gyration radius \/q cx iJ,~^^'^, we 
need /x cx N~^/'^. Thus, for a long polymer fj, ^ 0. We 
will also take /i = /io to simplify the algebra a bit. 

The pair of equations (85) and (86) fully determine jl 
as function of T and Tq. For /xq = M they take the form 



a(A) = 



+ 



1 



2^ 

1 



= + 



il{a) =fj,+ 



8TV/^(a + M-V2)f+i 



TTo{a 



a 2 



-1-1 



(87) 
(88) 



Given /i, T, Tq one can find the overlap with the native 
state <p and the size of the polymer from (81). 



A. Random-globule state 

It is immediately evident that when both the tem- 
perature T and the selection temperature Tq are large, 
/u « /U in (86), leading to a random-globule solution 
a = b = H^^^^, q = ^^^/^/2, ip = 0. What is not so 
obvious is that in the n ^ Q limit a solution very close 
to this exists all the way down to very low temperatures, 
even for small Tq. In this subsection we examine this 
state in detail. 



We look first for solutions of Eqs. (87) and (88) with 
the ansatz a = fi/fj, fixed and ^ 0. We call this the 
random globule ansatz, since, as will be shown, the poly- 
mer docs not have any fixed conformation (it is melted), 
and on the average the conformations it adopts have zero 
overlap with the native state. (Strictly speaking, this is 
the only truly ergodic phase we find.) For a we get, 



1 



3+i 



4^/5 



(89) 
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which, after inserting into (88), gives 



a : 



(90) 



Eqn. (90) can be used to calculate a as a function of /x. 
One can sec easily that a ^ I when /i 0. This shows 
that our ansatz is self-consistent in the limit of small fi. 
Also, (84) and (81) become 



1 



q = 



2^ 

1 



a 



1 



[l + 0(a-l)]. 



(91) 

(92) 



The normalized overlap between the polymer confor- 
mation and the native state is: 



COS0 



lim,_^(.r(^.O.r„(^)) 



Vlimt^oo(x(s,t)2)(a;o(s)2) ^q{^) 



■ (93) 



^From (91) and (92) we get cos(6») - (a - 1) - 
Thus there is no overlap with native state as — > 0. 

Furthermore, to check that polymer does not freeze 

into some other conformation, we calculate the normal- 
ized overlap between two configurations taken at very 
different times, 



cost^ = = = -. 

V[linit^oo(a;(s,t)2)]2 q 



After rewriting 



q/q = 1 ^ = 1 ^= , 

^'^ 2q 2V^' 



(94) 



(95) 



4TTn 



1 < 1 



(96) 



for the existence of a random-globule-like state. 

Some caution is in order. Working this out for finite 
N, d = 3, and an average density of 1, we find that the 
inequality (96) is violated below a temperature 



7r\i/2 15^2 



8To 



Ar-i/3. 



(97) 



With the small power of A^"-'^, one has to go to quite 
large A^ to make this temperature very low. Thus our 
statement that the random-globule-like state exists for 
all temperatures in the fi ^ limit may be of limited 
relevance for real 3-dimensional heteropolymers of the 
length of typical proteins. Nevertheless, here we are just 
considering this simple limit. 

We now discuss the stability of this solution. In the 
large- A' limit, it is locally stable against spontaneous for- 
mation of a native- like state at any T and Tq. However, it 
is unstable against glass formation at low temperatures: 
Since it is identical with the random-globule solution of 
the completely random hetcropolymcr problem, we can 
take over the result from paper I that it is unstable below 
a temperature Tg oc B, with the constant of proportion- 
ality of order 1. This glass temperature is independent 
of Tq. (In Fig. 1 this is the transition at /i ^ 0.) Thus, 
wherever the system is in a random-globulc-like state at 
T > Tg, it will no longer equilibrate if the temperature 
is lowered below Tg. Instead, it will become glassy and 
its dynamics will show aging. 



B. Ergodic native state 



and, using (92), we get cos^' = 0{a — 1). Again, as 
H ^ 0, cos 9' 0. This confirms that the ansatz 
a = 0(1) and /i leads to a melted random-globule- 
like phase. This phase is identical to that found at high 
temperatures for the completely random heteropolymer 
in paper I. 

The validity of the present ansatz rests upon the fact 
that we can solve Eq. (90). Clearly, for /i ^ a solution 
can always be found, namely a = 1. Since the physi- 
cally relevant is ex N~^/'^, we can always satisfy this 
equation, for any Tq, in the limit N ^ oo. 

We now address briefly the question of what happens 
for finite A^ (and /i). One can easily sec that Eq. (90) has 
two solutions when /(TTq) is not too large (e.g. 

by plotting the left- and right-hand side as functions of 
a) . The solution close to 1 is lost when the slopes of the 
left- and right-hand sides become roughly equal. Evalu- 
ating these slopes leads to the condition 



At low To and T, one expects that the polymer should 
be very close to its native state, i.e., small a. Therefore 
we also look for such solutions of Eqs. (87) and (88). We 
will try to solve equations (87) and (88) in the limit where 
/U ^ and jl stays finite. The limit fj, ^ turns out not 
to involve any subtleties when jl is kept constant, so we 
will just set /i = from the outset. Eqs. (87) and (88) 
become 

3 

""^^^ " ^ 8TV/'(a + /i-V2)*+i ^^^^ 

m = (99) 

Tro(cr-t-a)2+i 

These equations can be solved for fi as function of T and 
Tq. However, one has to keep in mind that n ^ has 
been taken. This implies that (81) and (84) become 
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and, inserting (100) into (93), the normahzed overlap be- 
tween native state and polymer conformations, becomes 
cos 9 ^ 1. Furthermore, because of its large overlap with 
the native state, the polymer is essentially frozen. This 
can be seen by calculating the normalized overlap be- 
tween two polymer conformations after a very long time 
interval, as in previous section. Inserting (100) into (94) 
and (95) gives cos^' Ri 1 — ~^ 1- 

There is interesting behavior associated with the limit 

fi ^ for very long polymers. When the polymer gets 
longer and longer [N — !■ oo) a finite part of the chain is 
not in the native state conformation, since a stays con- 
stant. The rest of the chain is in the native state, which 
can be seen from the fact that overlap with native state 
approaches 1. Thus, in the limit of a very long polymer, 
the fraction of chain not in the native state conformation 
becomes negligible: the recipe for biasing the coupling 
constants B^s' described in chapter II works best for long 
polymers. 

In the following we will proceed with the solution of 
equations (98) and (99). Before continuing, it will be 
useful to cornpactify notation a bit. Making the change 
of variables X = X/a for X = h^, a, q, q;Y = Ya^ for 
Y = H,^; and Z = Za^'^'^^/^/B for Z = T,To,we get 
equations of the same form, with X ^ X , Y ^ Y and 
Z ^ Z, but with (7 = 1 and B = 1. Thus, without loss 
of generality, we can choose units with a = 1 and B = 1 
(and remove the hats). From now on we do this. 

The working strategy for solving the equations is as 
follows. For fixed T, one can consider To as a function 

of p,. This can be easily done by inserting the expression 
for a from Eq. (98) into (99), thus writing /3o = 1/Tb as 

/3o{^T) = Tjl[l+a{il,T)f^+' (101) 

The four panels of Fig. 2 shows the shape of Po{jj,,T) as 
a function of fl for four different temperatures. We want 
ultimately to construct a phase diagram in the {Po,T) 
plane. Therefore we have to specify T (one panel of the 
figure) and /3o and ask whether one or more solutions, 
i.e., particular values of jl which solve Eq. (101), exist. 
For example, in panel (a) in Fig. 2, a horizontal line at 
00 > /3(f'" intersects l3o{ix,T) curve at two places, in- 
dicating two solutions jl = fii, jl2- To make the figures 
more readable we have shown such a horizontal line, at a 
particular values of /3o, only in panel (a). If this horizon- 
tal line is moved below /3™'", it will never intersect the 
/3o(/u,T) curve. Thus, we can see that for every T, there 
is a value /3™'"(T) below which no solutions exist. 

We proceed with the analysis of Fig. 2. For sufficiently 

high temperatures (panels (a)-(c)) there are exactly two 
solutions for all Po > /3™'". Of these, the one with the 
larger value of /i is a stable solution (local free energy 



minimum) describing the ergodic native phase. For ex- 
ample, the solution labeled /i2 in panel (a) is of this sort. 
The one with the smaller value of fl (e.g., the one labeled 
by /Ui in panel (a)) is unstable. It describes a free energy 
maximum between the minima at the random-globule 
and ergodic native states. We will call such states "un- 
stable stationary" (abbreviated US). (We have not done 
a static calculation to show this, but the situation here 
is analogous to that in an ordinary ferromagnet below 
the Curie temperature. There, one has three solutions 
of the mean field equations, one with positive, one with 
negative, and one with zero magnetization. The middle 
one, with zero magnetization, is unstable). The US state 
has a lower overlap with the native conformation than 
the ergodic-native solution does, because it has a smaller 
value of jl. As (3o is increased from below through 
the native-state and US-state solutions appear together 
and separate. For the temperatures of panels (a)-(c), 
they both exist for all /3o > /SJT'''- 

Panel (d) (at the lowest of the temperatures) shows a 
more complex behavior where double- minimum structure 
appears. We have found numerically that this happens 
below T w 0.20. Here the behavior around is just as 
in the other cases, but we note that at this temperature 
(3o{il,T) has a second local minima at a smaller value of 
fl. Thus there is a range (i™™ < Pa < Z?™'"'' for which 
there are four solutions. The rightmost one is stable and 
describes the ergodic-native phase, as before. Moving 
from right to left, the solutions alternate between stabil- 
ity and instability. Thus the second solution from the left 
represents a locally stable conformation. It is also corre- 
lated with the native state, since jl is finite (though we 
always find /i <C 1 in 3 dimensions). The remaining two 
solutions (with dPo/djl < 0) represent US states (local 
free energy maxima) between it and the random-globule 
phase in one direction and the ergodic-native phase in 
the other. 

Plotting against T, we obtain the stability bound- 
ary indicated by the thick solid curve in Fig. 3. Within 
our present assumption of ergodicity, everywhere to the 
right of this line the ergodic-native phase is dynamically 
stable. One can invert the relation /3™'"(r), obtaining 
a transition temperature T„(/3o), the maximum temper- 
ature for which the ergodic native phase is dynamically 
stable. It is separated from the (also stable) random- 
globule phase by a barrier, the top of which is described 
by the unstable solution. 

In Fig. 3 we also indicate the region in the {Po,T) 
plane where the second locally-stable solution is found. 

This region has the form of a kind of sliver extending out 
toward large /3o at low temperatures. 

So far we have not examined the stability of these so- 
lutions against glassiness. As indicated above, we do this 
with the help of Fig. 1: Stable solutions can not lie in 
the range flmin < A < P-at- In Fig. 2, these limits are 
marked on the jl axes. We thus see, for example, in Panel 
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(c) , that the native-state solutions found for the range of 
Po corresponding to values of jl between /i* and fiAT are 
not acceptable: they violate the AT stability condition 
(82). 

Similarly, in panel (b) the US solutions found for a 
range of f3o values can also be seen to lie in the forbidden 
region. And the intermediate locally-stable states that 
we identified in panel (d) always lie in a glassy region. 

In Fig. 3 we also plot the AT Une (80) in the {/3o,T) 
plane, indicating the regions where the various kinds of 
ergodic solutions are forbidden. For the native-phase so- 
lutions, the forbidden region is a strip mostly at low val- 
ues of T (diagonally cross-hatched region between thick 
and AT line). However, it "wraps around" at the leftmost 
part of the region where those solutions are found. 

The forbidden region for the US solutions occupies 
most of the region where these solutions occm below 
Tmaxi the maximum temperature for a glass transition 
shown in Fig. 1, including the entire portion of it be- 
low Tg, the glass instability temperature of the random- 
globule state. 

The structure in a tiny region near the minimum value 
of /3o for which ergodic-nativc solution arc found is a bit 
complicated and cannot be seen in the top panel of Fig. 3. 
Therefore, the lower panel shows an enlargement of this 
region. 

In summary, we have found four kinds of ergodic solu- 
tions. One essentially describes a random globule state. 
It is locally stable (in the limit of a large globule) at all 
T and (3q against condensation into a native-like state, 
but unstable against glass formation everywhere below 
a transition temperature Tg. The second kind of solu- 
tion describes a phase which is highly correlated with 
the native state conformation, and it is stable in most of 
the region where the solution exists. The third kind of 
solution describes a locally stable state, correlated with 
the native state but more weakly so than the ergodic na- 
tive phase just described. It is never stable against glass 
formation. Finally, there are unstable solutions, found 
whenever the ergodic-native solutions exist. They de- 
scribe US states, free energy maxima between pairs of the 
previously-described solutions. However, in a large part 
of the region where these solutions are found (roughly, 
everywhere below T^ax ~ Tg) they violate the AT sta- 
bility condition and so are not physically relevant. 

Outside the regions where these ergodic solutions are 
allowed, we have to look for glassy solutions. We do this 
in the next section. 



X. GLASSY PHASES 



In a glassy phase, f(l) ^ and Eq. (80) has to be kept, 
which gives, 



= -6'^V"(6) 



(102) 



Also, equations (74), (75) and (76) can be rewritten in 
the form 



V{h)-V{ho) _T^^(l 



bo-b 



2 b \b 



bo-b=- i-^-b 



and, with fiQ = fi, (78) and (79) become 

_ b_a 1 J_ 

" ~ 2 ^ 8TV/^(1 + bo)i+^ ^ 
1 



(103) 
(104) 



1 



TTo(l + a)i+i 



(105) 
(106) 



The above equations can be solved as follows. Eq. (102) 
gives & as a function of T, and then (103), (105) and (106) 
can be used to find 6o and Jl as functions of T and Tq. 
Once 6o and jl are found one can calculate q as 



1 



2 8rv/2(i + j,g)i+i 

4^' Vm 



(107) 



As in our analysis of ergodic solutions in the preced- 
ing section, we will try two types of ansatz: one with 
jl/ 11 = const as ^ — > and one with jl = const as /U ^ 
leading to what we call frozen-globule and glassy native 
phases, respectively. 



A. Prozen-globule phase 

The limit where a = /i/yit is kept constant and /U ^ is 
easily treated. Eq. (102) stay the same, while Eq. (103) 
gives 



V{b) - V{bo) T^yS 



bo-b 



262 



(108) 



Since b is kept fixed the only solution of equation above 
is 6o ^ oo as 



bo 



afj, 



(109) 



where V(&) is a function which depends only on 6, as fi 
is sent to 0. Inserting (109) into (106) gives 



a = 1 + 0{ti^''-^y^) 



(110) 
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and a stays very close to 1, as in the ergodic random 
globule case. Also, is given by (91), while (107) gives 



1 



2VM 



(111) 



which can be compared with ergodic globule result, 
Eq. (92). Eq. (93) stays the same, and one gets cos^ 
a — 1 which goes to zero as /x ^ 0. There is no over- 
lap with native state. Does the system freeze into some 
other configuration? To find out, we calculate overlap 
angles between configurations at time t and a much later 
time t' . As discussed in section VII there are two ways in 
which the limit t — > oo can be taken, leading to qo ^ q. 

In the first limit, the equivalent of Eq. (94) for the 
ansatz used here reads 



, \\mx,t^ao{x{s,t)x{s,\t)) qo 
cos6'„ = =^ = — 



(112) 



: 1 — &o/2(7, and 
1 which goes to 



Using Eq. (64), we can write cos^g 
Eqs (111) and (109) give cos6»^ ~ a 
as /i — > 0. (This behavior is analogous to that found 
in p-spin glasses.) However, at not-too-long time scales 
(shorter than the waiting time), as in equation (113), the 
polymer is frozen: 



cos6lg 



lim. 



■T,t—^ OO 



{x{s,t)x{s,t + t)) 



V[limt^oo(a;(s,i)2)]2 



(113) 



Using Eq. (64), we can write cos^^' = 1 — b/2q, and 
Eq. (Ill) gives cos^g 1 — ^b/iplb), which goes to 1 
as /X ^ 0. Thus, this glassy phase has no overlap with 
the native state. 

As discussed above, there is an upper temperature 
limit Tg (independent of /3o) above which this phase 
melts, leaving the system in the random- globule state. 
Tg can be found from Eqs. (102) and (103), using 
&o ^ oo and (104) with x — > 1. This leads to a 
value b — 2/{^d — 1) = 0(1) at the transition and 

Tg = 2(id - l)^(l'^-iV(5d + l)5(5'^+i). For d = 3, 
T„ w 0.535. 



B. Glassy native states 



We also have to study the possible glassy phases with 
overlap with the native state, i.e., with finite p, (and, ac- 
cordingly, finite a) when /x ^ 0. In such a phase, as in 
the ergodic native-like states described above, the sys- 
tem moves only in the neighborhood of the native state 
configuration. However, in a "glassy native" state even 
this restricted motion is strongly suppressed by the com- 
plexity of the local potential energy surface, and a glassy 
phase results. 



As in the ergodic ansatz, the limit /i ^ introduces 
no problems. Eqns. (102) and (103) remain the same as 
in the frozen-globule case, while the equations for a and 
jl, (105) and (106) become 



a(M) = TT + 



fl{a) = 



1 



8^2^3/2(1 
1 



bo) 



+1 



4Vm 



TTo(l + a)S+i 



(114) 
(115) 



Again, Eq. (102) specifies 6 as a function of T, and (103), 
(114) and (115) determine bo and jl as functions of T and 
Tq. q and ip are given hy (p,q 1/(2^//!). 

The overlap with the native state is the largest possi- 
ble: cos^g = 1, as can be easily seen from Eq. (93) and 
the values for ip and q we have just given. The over- 
lap between two conformations at very different times 
also takes its largest possible value. From Eqs. (112) 
and (113), knowing that bo and b do not depend on 
/X we have cos^g = 1 — bo/2q « 1 — bo^Jl 1 and 
cos6'g = 1 — 6/2(7 « 1 — by^ 1. Thus, the polymer is 
frozen almost everywhere into the native conformation. 
However, the freezing is not total, since a in (114) is not 
zero. Furthermore, there is aging in the system, since x 
in Eq. (104) is not equal to 1. 

We turn now to the solution of the equations (102), 
(103), (114) and (115). As for the corresponding ergodic 
phases we have to resort to numerical solution; here we 
describe the analysis. The working strategy is similar to 
the one presented in subsection IX. B; the goal is to find 
/3o as function of fi for fixed T since, as in the ergodic 
native case, extrema of the function l3o{jl,T) govern the 
phase boundaries. 

The procedure for finding value of the function 
To{jj,,T) is as follows. Eq (102) determines 6 as a func- 
tion of T, to be referred to as b{T). Once b(T) is found 
from (102) it is inserted into Eq. (103), which deter- 
mines 5o(T, /i). The value found for bo is inserted into 
Eq. (114) to find a, and finally /3o = l/7o is calculated 
from Eq. (115). Thus, at each temperature for which 
glassy solutions are possible, we can construct a graph of 
/3o(/i); as we did for ergodic solutions in Fig. 2. We have 
used Mathematica to do these calculations. We can use 
these curves, together with the ergodic ones previously 
analyzed, to identify the possible states of the system at 
a given temperature and Po (Fig- 4). The procedure is 
fairly simple. At any given jl, only one of the solutions is 
physical: In the region jlmin < jl < iIat one has to follow 
the glassy /3o(A) curve, while outside it one follows the 
ergodic one. In Fig. 4 the physical solution is indicated 
as the thick dashed curve. One then looks for solutions 
as intersections of this curve with a horizontal line at a 
given value of /3o, as done previously (e.g., as in Fig. 2, 
panel (a)) within the ergodic ansatz. 

In Fig. 4 this procedure is shown for several different 
values of T. In the first panel, T lies just a little below 
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Tmax (as in Panel (b) of Fig. 2). Suppose we start in 
the ergodic native phase at largo f3o and then lower Po- 
(In Fig. 5, this would correspond to moving along a hor- 
izontal line (constant T) slightly above point B in panel 
(b) or (c)). We can lower f3o all the way down to 
without cnc;ountcring an AT instability. So, just as in the 
ergodic analysis of section X.B, beyond the ergodic 
native phase melts into the random-globule phase. 

In the same panel we can also analyze the what hap- 
pens to the unstable stationary state in the same range 
of /3o for this temperature. At very large f3o we have an 
ergodic solution, but as we lower f3o we pass through a 
range of //, between flmin and jlAT, where the ergodic 
solution is unstable against glassiness. In this region 
we must follow the glassy curve instead of the ergodic 
one. We interpret this glassy solution in the following 
way: The free energy landscape near the US maximum 
becomes rough in this range of values of /3o (at this tem- 
perature), the same way the free energy landscape near 
the minimum corresponding to a thermodynamic phase 
becomes rough in a glassy state. We call it a "glassy US 
state" . 

The next panel is for a slightly lower temperature (but 
still above Tg). Here, as we lower f3o in the ergodic- 
native phase, we reach an AT instability before we get 
all the way down to the minimum on the ergodic curve. 
(In Fig. 5, this would correspond to moving on a line 
of constant T, meeting the AT line somewhere between 
points A and B in panel (b) or (c)). Furthermore, the 
only available glassy solution for fl < (Iat is one with 
negative d(3o/dfi, that is, it corresponds to the kind of 
glassy US state discussed above. As this is not a stable 
phase, we conclude that for this T, the minimum value 
of /3o lies at this AT line, and beyond it there is no stable 
native-like state. We can follow the glassy US state back 
up to larger /3o, seeing that we eventually cross over to a 
normal (non-glassy) transition state. 

In the last panel, the temperature is lowered a bit more 
(below Tg). Again, starting in the ergodic native phase 
at large /3o and lowering /3o, we encounter an AT insta- 
bility and a glassy solution appears. (Equivalently, in 
Fig. 5 one moves on a horizontal line somewhere below 
Tg until meeting the AT line for the first time.) For 
smaller /3o, we switch to the glassy curve, which has pos- 
itive dPo/djl, describing a glassy native phase. We can 
follow this curve down to its minimum /3o, beyond which 
no phases correlated with the native phase exist. But, 
of course, following it back up toward large /3o on the 
unstable branch, we can identify the glassy US state be- 
tween the phase correlated with the native state and the 
one uncorrelated with it. (Above Tg, the latter is the 
random-globule phase; below it, it is the frozen-globule 
phase.) 



XI. FEATURES OF THE PHASE DIAGRAM 



The phase structure implied by this simple model is 
not so simple. Fig. 5 shows the phase diagram con- 
structed from the above analysis. For clarity, we show 
in the top panel only the solutions that correspond to 
stable phases. The second panel shows the details in the 
region where the ergodic-native, glassy-native, random- 
globule and frozen globule states come together (or nearly 
so) . The third panel shows the regions where ergodic and 
glassy US states are found. 

There are six distinct regions in the phase diagram. In 
region I (high T, small /3o), the only stable phase is the 
random globule. In region II (small 0o, T <Tg) it under- 
goes a glass transition to the frozen-globule phase. The 
properties of the system in this part of the phase diagram 
are the same as in the completely random heteropolymer 
model of paper I; the bias of the interactions toward a na- 
tive state does not have any effect until a (temperature- 
dependent) threshold PgiT) is reached. This threshold is 
marked on the diagram by the lines separating region I 
from regions III and V and region II from region VI. 

To help thinking about these phases, we offer the 
schematic free energy-surface pictures of Fig. 6. They 
show how we imagine the free energy varies as a function 
of the native-state overlap coordinate (p. Fig. 6A depicts 
this cross-section through the free-energy surface in re- 
gion I, where there is a single smooth minimum around 
ip = 0, representing the random- globule phase. Fig. 6B 
shows what happens in region II, where this phase is re- 
placed by the frozen-globule phase. We represent this 
by drawing the free energy surface with many local min- 
ima. Fig. 6C shows what happens in the middle of region 
III, where there are two (smooth) minima, the new one 
corresponding to the ergodic native phase. (It will lie 
above or below that at tp = according to whether we 
are above or below a first-order transition that we ex- 
pect to occur at a temperature Ti < T„, see below.) In 
Fig. 6D we depict the situation in region IV, where both 
the Ri region and the region aromid the maximum 
become rough. Fig. 6E represents region V, with the na- 
tive valley and the maximum rough, but the region near 
(f = still smooth. In Fig. 6F (region VI), that, too, 
becomes rough. 

An important feature is the fact that the random- 
globule and frozen-globule phases remain (dynamically) 
stable in their respective temperature ranges for all /3o- 
Thus the horizontal line separating region I from region 
II continues across the diagram, separating region V from 
region VI and region II from region IV. 

In each of regions II, IV, V and VI (Fig. 6, panels (b), 
(d), (e), (f)) there are two stable states. Thus, regions I 
and III arc separated by the line r„(/3o)j below or to the 
right of which, in addition to the random-globule state, 
an ergodic native state exists and is stable. In region 
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IV, the frozcn-globulc and this crgodic native state arc 
both stable. In region V, the random globule state and 
a glassy native state are stable, while in region VI the 
stable states are the frozen globule and the glassy native 
phase. 

The diagram shows some interesting fine structure in 
the neighborhood of the region where regions I, III, and V 
meet (second panel, point A). In particular, below point 
A, the boundary between regions III and V (i.e., between 
the ergodic-native and glassy-native state) is an AT line. 
It comes about as can be seen in the last panel of Fig. 4, 
where, following the ergodic phase down from high (3o 
and fi, we reach jiAT and thereafter have to switch to the 
glassy solution. 

On the other hand, the boundary r„(/3o) above point 
B is reached as in the first panel of Fig. 4: one can come 
all the way down to the minimum value of /3o for the 
ergodic solution before reaching fiAT- The full line con- 
tinuing upwards and to the right of point B is an AT line 
which goes over into an a; = 1 line at its maximum, T^ax- 
Below this line US solutions become glassy. 

Between points A and B, the boundary is marked by 
reaching JIat in the way shown in the second panel of 
Fig. 4: There is no stable ergodic native solution to the 
left of this line, since, upon lowering /3o below in 
Panel (b), the solution with d(3o/dfi > is lost. Thus in 
this region the line AB is an AT line for both the native 
phase and the US solutions. The dashed line marks the 
boundary found if one ignores the AT line (i.e., it is a 
portion of the boundary found using the ergodic ansatz 
and shown in Fig. 3). 

Ever}rwhere below the AT line (both the portion AB 
and its extension upward to the right) is the region where 
the US solutions are glassy, as shown in the last panel of 
Fig. 5. At a given /3o, these features all have an onset 
at a temperature between Tg and Tmax- As can be seen 
in Fig. 1, this is a very small temperature range. This is 
also the reason why there is fine structure, as shown in 
both the second and the third panels of Fig. 5, in such a 
small temperature range in the phase diagrams. 

As remarked above, we have not done an equilibrium 
calculation, but we expect that the first-order transi- 
tion temperature Ti{f3o) where the free energies of the 
random-globule and ergodic-native phases are equal will 
also rise with /3o- For large /3o, we expect Ti, like T„, to 
be proportional to /3o (but Ti < T„, of course). 

Thus, at fairly large (3o, we expect the following se- 
quence of stable states as we lower T from a high value. 
Initially, only the random-globule state is stable. Then, 
below Tn, the ergodic native state is also stable, and be- 
low Ti it becomes the lowest-free energy state. Going 
further down in T, we cross the boundary (last panel in 
Fig. 5) where the US state between the ergodic native 
and random-globule states becomes glassy (i.e., acquires 
a rough local free energy landscape). Very soon there- 



after, we cross Tg, where the random-globule state under- 
goes glassy freezing. Continuing, we reach a temperature 
where the ergodic native state undergoes glassy freezing. 
Finally, we reach the stability limit of this glassy-native 
phase, leaving the system with nowhere to go but the 
frozen-globule state. 

What lessons are there in these findings for protein 
folding? We start from the assumption that the initial 
state in the folding process is uncorrelated with the na- 
tive state (i.e., in Fig. 6 we start in a local minimum at 
(fi 0). Folding requires the system to find its way to 
the (ergodic) native state. 

One feature that is evident is that such a path in 
configuration space always requires an uphill free-energy 
step. This is because either the random-globule or frozen- 
globule phases is always locally stable. 

If we stick to our mean-field dynamical picture, where 
barriers are infinite, folding is, strictly speaking, impossi- 
ble. In dynamical terms the "infinite" barriers translate 
into the fact that the equations which govern the motion 
of order parameters, tp included, have basins of attraction 
corresponding to the plots shown in Fig. 6. For example, 
starting from (p somewhere close to and given a free 
energy profile like that in Fig. 6C, dynamical equations 
will never carry ip to the large value describing the native 
state. On the contrary, (p will approach as time goes 
on. 

But, if we relax this assumption and imagine finite bar- 
riers (associated with local nucleation of a native phase 
[44- 46]), we may ask (informally) when activated motion 
over the barrier to nucleation is least hindered. We argue 
that the free energy landscape features present globally in 
our mean-field picture will also be relevant locally: when 
our calculation here finds a glassy US state, we expect 
that free energy surface near the true transition state will 
also be rough. Thus, from the preceding description of 
the phases and the transition states between them, we 
can see that folding should be easiest for large (3o in a 
window between Ti(/3o) and the upper boundary of the 
region where the US states become glassy (and passage 
across the transition region is kinetically impeded by the 
tortuous nature of the local free energy landscape) . The 
latter boundary lies, in turn, just barely above Tg, where 
the landscape in the (large) portion of the configuration 
space uncorrelated with the native state also becomes 
rough, further impeding escape from it. At still lower 
temperatures, things become even worse, first with the 
onset of glassiness in the native-like region of configu- 
ration space itself and finally with the disappearance of 
native-like solutions. But these features probably have 
minor consequences, since folding will already have been 
so strongly impeded by the effects (with onset near Tg) 
that tend to confine it in a region of configuration space 
uncorrelated with the native state. 
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XII. DISCUSSION 



of these: 



We have introduced what we might caU a generic model 
for a protein, based on what seems to us to be the sim- 
plest way to incorporate a tendency to form a native state 
in an otherwise random heteropolymer model. To make 
it possible to calculate typical properties, we follow previ- 
ous authors [2,24-26] and do not specify a particular na- 
tive state, but rather an ensemble of them, constrained 
only by chain-entropic constraints and confinement to 
the appropriate volume. This ensemble is characterized 
by the selection temperature Tq. Our model differs from 
previous ones in that they are based on random-sequence 
heteropolymers, while we start from a model [6,7] in 
which each monomer-monomer interaction is an indepen- 
dent random variable. 

While it might be argued that random-sequence mod- 
els are more relevant to proteins, they approach the 
model we consider here in the limit where the number of 
monomer types becomes large. Thus, what we find out 
about our model may be relevant to proteins (with 20 
different amino acids). Of course, it is also important to 
study what happens away from the large-monomer-type 
limit; our analysis here can help in solving that more 
difficult problem. 

Furthermore, naively, one might assume that by ad- 
justing N(N — l)/2 parameters one could imprint a na- 
tive state more strongly than for models with only N 
parameters. Our model shows that this is not necessar- 
ily true. Parts of the phase diagram are glassy, even for 
very low selection temperature Tq, when the native state 
should be strongly imprinted into the model. 

Instead of the quadratic confinement term j2x{s, t)^ one 
could add three-body terms, which are commonly used to 
fix the globule density. It would be interesting to extend 
the analysis presented here to such models. Also, in our 
treatment, translational invariance within the globule is 
put in by hand. Keeping three-body terms would lead to 
automatic translational invariance. We have seen that if 
translational averaging is omitted (see paper 1) then the 
equations become coupled in the k variable and are thus 
a lot harder to solve. 

Within out model, we have made just two approxima- 
tions: the Gaussian variational ansatz of section VI, and 
the assumption of 1-step ergodicity breaking (analogous 
to 1-step replica symmetry breaking in the replica ap- 
proach). Otherwise the solution is complete and exact to 
the accuracy we were able to achieve numerically. 

Our most important result is the existence of the vari- 
ous different phases at large B /Tq, where the interactions 
that favor a native state are strong. While it is natural to 
anticipate that the native-like configurations will be ther- 
mally disrupted above a temperature of order B'^/Tq, it 
is not so obvious that at low temperatures there will be 
other impediments to efficient folding. We identify two 



(1) The frozen-globule state, which is uncorrelated 
with the native state, always exists below Tg, no matter 
how big po is. This means that in a large part of con- 
figuration space, the system may be trapped in a rough 
energy landscape and never (in MFT) get to the native- 
state region where it can fold rapidly. Furthermore, in 
almost the same temperature, range, we expect that the 
energy landscape is also rough around the transition re- 
gion on the way to the correctly-folded state, further im- 
peding the folding process. Thus, while lower tempera- 
ture favors well-folded over random-globule-like configu- 
rations energetically, the rough energy landscape of the 
glassy phase will hinder correct folding. Our conclusion 
here is consistent with that of Goldstein et at. [47], who 
found (albeit in a different kind of model) that a large 
Tn/Tg (or Ti/Tg) ratio favors folding. 

(2) At even lower temperatures, the native state itself 

is unstable against a glass transition where it splits into 
a large number of substates. Transitions between these 
substates are blocked by high barriers (infinite, in MFT). 
A phase of this kind was found earlier by Bryngelson and 
Wolynes in a phenomenological model [48] . ft is tempting 
to associate the substates with the glassy conformations 
observed at temperatures below 200 K in myoglobin [49] . 

Of course, MFT is an approximation. The escape from 
the tortuous part of the energy landscape to the smooth 
region will not take forever, nor will transitions between 
low-T substates. Nevertheless, MFT does indicate when 
we can expect relaxational dynamics, including folding, 
to be slow or fast, as well as give us some insight into the 
physics behind these differences. 

Our analysis here is a purely dynamical one. We do 
not compute equilibrium partition functions. A com- 
plete analysis would include such calculations, but we 
defer them to future work. Nevertheless, the purely dy- 
namical analysis can reveal important properties of the 
system that cannot be seen in an equilibrium analysis. 
For example, it has been know for a long time that for a 
large class of models — namely, those which have a glass 
transition where x ^ 1 — the dynamic and static glass 
transition temperatures are different [39,42,43]. This is 
expected to be the case for the transition at Tg in our 
model: The equilibrium glass transition temperature is 
lower than the dynamical one. Thus, in a temperature 
range just below the dynamical Tg, the equilibrium anal- 
ysis does not reveal the slow dynamics (accompanied by 
aging) that we are able to identify and analyze here. 

In other glassy models for which it has been possible to 
do a more complete analysis [34,35,39,42,43], the static 
and dynamic transitions coincide when they occur as a 
result of an Almeida-Thouless instability (the marginal 
stability condition, Eqn 80) [50] . This is the case here at 
the phase boundary where the native-like state becomes 
glassy. 
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Gillin and Sherrington [51] and Gillin et al. [52] have 
been able to analyze both the statics and the dynam- 
ics of several classes of mean-field spin glass models with 
a competition between glassy and ferromagnetic states 
(see also ref. [53] for a special case). Some features of the 
phase diagram of our model that we have been able to 
discover so far are also seen in their models. 

Gillin et al. also studied full (as well as 1-step) replica 
symmetry breaking solutions, which we have not. In 
some of their models, the counterpart of our glassy native 
phase undergoes full RSB at low temperatures, and the 
counterpart of our II- VI boundary becomes vertical. It is 
possible in our model as well that, in particular regions 
of the phase diagram, our 1-step solutions are not sta- 
ble and full RSB is necessary. More genc;rally, it will be 
an interesting problem to try to explore what kinds and 
degrees of universality there are in the phase structures 
of various systems where the glassiness induced by disor- 
der competes with some kind of order analogous to the 
native state in our problem or the ferromagnetic state in 
theirs. 



$. (The same reasoning holds for a transformation like 
(A2) with a derivative with respect to instead of with 
respect to 9.) The only term left is the S'i[<i>] (i.e., the 
part of the action quadratic in the superfield) and it is 
straightforward to see that this term is also invariant un- 
der (A2) (though not under the transformation involving 
the derivative with respect to 9). 

The fact that the action is invariant under (A2) implies 
the Ward identity 



d 



on. =0 



(A4) 



which gives = {v^o) + 9i{xxo) (we have suppressed the 
arguments of the fields to simplify the notation). This 
implies that separately one has 



{rixo) = 0, {xxo) = 



(A5) 



Inserting (A5) into (Al) gives the desired result, Eq. (30). 
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APPENDIX A: CORRELATION FUNCTION Gff 

Here we derive Eq. (30). Inserting Eq. (14) for the 
superfield into (25) gives 

GfS' = {x{s,h)xo{s')) + {f]{s,h)xo{s'))ei + 

+ei{v{s,h)xo{s')) + eiei{x{s,h)xo{s')). (ai) 



Here we give more background on the use of Eq. (36). 
In the dynamical calculation, the fields are complex and 
contain Grassmann variables; thus. F^yn is not a real 
number. This means that any interpretation of Eq. (36) 
as an extremum condition for Fdyn has to be given up. 
Nevertheless, we can still make some sense of the GVA 
as the first step in a systematic approximation scheme. 

Formally, one starts from Eq. (23), which we rewrite 
in the shorter form 



(Bl) 



One can show that the action of the dynamical generat- 
ing functional Fdyn (see Eq. 23) is invariant under the 
infinitesimal transformation (BRS symmetry) 



(A2) 



This follows in two steps. First one notices that for any 
function / 

S I Sfim) = e J Stfi^e)) = 0, (A3) 

due to the identity J d9-^ = 0. This means that any 
term involving a local function of $ (i.e., not containing 
derivatives over 9 and 9), e.g., S2[^,xq] (see Eq. 16), is 
invariant under the transformation (A2). S'o[a;o] is triv- 
ially invariant since it does not contain the superfield 



where ^ stands for the pair {xq,^) and, likewise, D'^ 
for DxqD^. One can express F^yn in a slightly different 
form 



where 



and 



Fdyn = ^g-(S-S„ar.) 



fvar^ 



((•••)>^ 



/jjvl/(...)e-^- 



(B2) 



(B3) 



(B4) 



In a static calculation one would proceed with the in- 
equality 
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(B5) [s; 



to conclude that 

F ^ (^S Suar) var H" -^t;ar- 



(B6) 



Thus, in a static calculation, the variationally-obtaincd 
F gives an upper bound on the true F. What is allowed 
to vary is the form of S^ar, most often, the parameters 
describing it. (In the GVA, Svar is specified by Gfl and 

In the dynamical problem we follow another route, 
starting exactly at the problematic step, Eq. (B5), along 
the lines of ref. [54]. Instead of the inequality (B6) we 
use Eq. (B2) in a slightly modified form 

^dyn ~ ^var In^e }varj (■^'^) 

where AS = S — Syar- Applying a cumulant expansion 



(exp(-AS'))„„r= exp [-{AS)^ar+ 



(AS) 



var J 



one gets 



F'dyn — -^var ~l~ ^var) var ~l~ 



, (B8) 



(B9) 



where AF contains second and higher order corrections 
in AS. In any approximation made by keeping a fi- 
nite number of terms in (B8) (the simplest being to set 
AF = 0), F^yn depends on Gfl . To minimize this de- 
pendence, we chose G\2 so that the derivative of the 
approximate form for Fdyn with respect to Gfl vanishes. 
This gives Eq. (36). Furthermore, if all terms in AF are 
kept, this procedure, by construction, formally gives back 
the exact -Fdy„. 

The meaning of minimizing the dependence with re- 
spect to quantities involving Grassmann numbers may 
seem obscure, but we note that we arc using the SUSY 
representation only for compactness. The entire GVA 
calculation could have been presented equivalently with- 
out any Grassmann variables, with no change in mean- 
ing or result. Thus, we are really only minimizing the 
dependence on parameters of physically well-defined cor- 
relation and response functions. 
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FIG. 1. Boundary of the glassy phase in the {p.,T) 
plane. The boundary is same as in the case of the ran- 
dom-heteropolymer model from paper I, except that na- 
tive-state correlations lead to the replacement of /it by jl. We 
have used parameters d = 3 and a = 1. Tmax is the max- 
imum T for which Eq. (80) has a solution (see Paper I for 
further details). Hc is the value of ft, where T{jl) attains this 
maximum. The solid part of the boundary is an AT line, and 
the dash-dotted part marks a transition where a; — > 1. Ap- 
proaching the AT line from below, 6 — bo 0, while x remains 
strictly less than 1. Approaching the x = 1 line from below, 
a; — > 1 smoothly, while 6 — 6o is discontinuous there. Above 
both lines, b = bo and x is undetermined (any x ^ solves 
(74)-(80), and no physical quantity depends on it). The same 
holds for all figures where these lines appear. 





jl .=0 

' min 



AT \^ 



FIG. 2. Analysis of ergodic-native solutions: l3o{jl,T) as 
function of jl for four values of T (see text for explanation). 
Panel (d) shows the existence of the two extra solutions (one 
stable, the other unstable) in the range [PT'^ , pT'"]. 
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FIG. 3. Regions of existence and stability of ergodic-native 
solutions. Panel (a): Ergodic-native-phase and US solutions 
exist everywhere to the right of the thick solid curve. The 
ergodic native phase is stable against glassiness everywhere 
there except in the diagonally cross-hatched region. The US 
states are also unstable against glassiness there, and addi- 
tionally in the horizontally cross-hatched region. The vertical 
cross-hatching marks the region where the extra phase seen 
in Panel (d) of Fig. 2 is found. (This phase is never stable 
against glassiness.) Panel (b): Enlargement of the circled re- 
gion in Panel (a) . The AT line is tangent to the ergodic phase 
boundary (thick line). At its maximum, at Tmax, it becomes 
an ,T = 1 line (dashed-dotted line, sec also Fig. 1). Lowering 
T from the wliite region into tfic fiorizoritally cross-hatched 
region results in two different types of transitions depending 
on whether one crosses the AT or the x = 1 line. In both 
cases the US state becomes glassy. 
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FIG. 4. Analysis of glassy-native solutions: f3o{jl,T) as a 
function of /2 for three fixed values of T (see text for expla- 
nation) . Full linos: I3q {jl, T) calculated within the crgodic 
ansatz (as in Fig. 2). Dashed lines: l3o{jl,T) calculated with 
the glassy ansatz. The actual curves vary with /i in a way that 
is difficult to plot in a useful way, so here we have distorted 
them in such a way as to make their qualitative form (num- 
ber and ordering of maxima and minima) evident. When the 
two curves cross at /i = fiAT, one has to change from the er- 
godic to the glassy solution (when approaching from p, = oo). 
Similarly, when p, = jlmin one has to go back to the ergodic 
solution. The thick dashed line indicates the physically rele- 
vant states (both stable phases and US states). 
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FIG. 5. The final phase diagram. Panel (a): Stable phases. 
Region I: Random globule is the only stable phase. Region 

II: Frozen globule is the only stable phase. Region III: Er- 
godic native and random globule phases stable. Region IV: 
Ergodic native and frozen globule phases stable. Region V 
(only visible in panel (b)): Glassy native and random globule 
stable. Region VI: Glassy native and frozen globule stable. 
The dashed line marks the boundary of the (unphysical) er- 
godic native state from Fig. 3, to emphasize that the phase 
boundary of the glassy native state (solid) does not coincide 
with it. Panel (b): Enlargement showing structure in the re- 



gion near T = T„ 



Tg and /3o = 1.45 (including region 



V). Below point B the boundary of region III is given by the 
AT line. Above point B the boundary is the ergodic-native 
stability limit (the uppermost line in Panel (a)). The continu- 
ation of the AT lie is shown as a dotted line (which turns into 
dash-dot x = 1 line). Panel (c): US states are ergodic in the 
vertically hatched region, glassy in the horizontally hatched 
region. The boundary above and to the right of point A is 
an AT line. Beyond the region shown the a; = 1 line falls 
off monotonically, and for /So ^ oo it approaches Tg. Be- 
low point A, the small-/3o boundary coincides with the line 
between regions II and VI in panels (a) and (b) . 









FIG. 6. Schematic free energy surfaces in different regions 
of the phase diagram. (See text for explanation.) (a): region 
I. (b): Region II. (c): Region III. (d): Region IV. (e): Region 
V. (f): Region VI. 
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